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1. Introduction 

>, 
^D , Solutions of the hydrodynamical Riemann problem were introduced to 

J^ ■ the numerical hydrodynamics by Godunov already in 1959 [Ij. After 50 

years their significance is hard to be overestimated. Especially in the nu- 
^__^ ^ merical relativistic hydrodynamics, most of the so-called high resolution 

;<— ^ ■ shock capturing schemes is based on the modifications of the original Go- 

f^ I dunov idea (for a good review on this issue see p]). 

In general, by a Riemann problem for a set of hyperbolic partial differen- 
tial equations we understand an initial value formulation, where the initial 
data consist of two constant states separated by a discontinuity in the form 
of a plane surface. An elementary introduction and some theorems on such 
5^ I a general case can be found in [3]. Since in his original paper Riemann was 

concerned with equations of motion of the inviscid fluid, it is also customary 
to refer to the Riemann problem as a problem in hydrodynamics f^. 

The special case of the relativistic Riemann problem, where the gas in 
both initial states is assumed to be at rest with respect to a chosen iner- 
tial frame, was first considered by Thompson in 0. The general Riemann 
problem with arbitrary velocities, but in one dimension only, was solved by 
Smoller and Temple [B] for the ultra-relativistic equation of state and by 
Marti and Miiller [7] for the perfect gas equation of state. The latter work 

(1) 



was generalized to the case of arbitrary initial velocities, also tangent to the 
initial discontinuity, by Pons, Marti and Miiller in [8j. 

At this point the relativistic hydrodynamics differ strongly from its New- 
tonian counterpart. The solution of the Newtonian Riemann problem is 
independent of the velocities tangent to the initial discontinuity and the 
whole problem can be treated in just one spatial dimension. In the rela- 
tivistic case all components of the velocity couple to other hydrodynamic 
quantities through Lorentz factors, and even the wave pattern of the solu- 
tion can depend on the tangent velocity. 

A solution analogous to that of Pons, Marti and Miiller has recently been 
obtained for ultra-relativistic equation of state by Pi§tka and the author in 
[9]. Due to the simplicity of the ultra-relativistic equation of state, this 
solution can be expressed almost entirely in analytical terms, and as such it 
can be used both to construct and test modern multidimensional numerical 
codes solving equations of hydrodynamics. 

Our original motivation for dealing with solutions of the Riemann prob- 
lem was to construct a numerical scheme that could be used in simulations 
of cosmological hydro dynamical perturbations in the radiation-dominated 
universe. It is however worth pointing out that the issue is not only numer- 
ical. Recently, Aloy and Rezzolla used a solution of the Riemann problem 
with non-vanishing tangential velocities discussed in \8\ to explain a purely 
hydrodynamical mechanism boosting astrophysical relativistic jets to large 
Lorentz factors |10] . 

The ultra-relativistic equation of state, commonly used in cosmology, 
has the form p = kp, where p is the pressure, p the energy density, and 
k E (0, 1) is a proportionality constant {Vk can be interpreted as the local 
sound velocity). On the other hand, most numerical codes are adjusted to 
equations of state depending explicitly on the baryonic (rest-mass) density 
n and the specific internal energy e (satisfying p = n + ne). This is, for 
instance, the case of the standard perfect gas equation of state of the form 
p = (7— l)ne, where 7 is a constant. Although the ultra-relativistic equation 
of state can be viewed as a limit of the perfect gas equation of state with 
n — )■ 0, e — >• 00 and ne = const, a numerical code has to be rewritten in 
order to evolve a fluid with the ultra-relativistic equation of state. 

In this paper we will summarize the analysis given in [9], and then de- 
scribe the construction of a 3-dimensional relativistic hydrodynamical code 
with the ultra-relativistic equation of state. We will also show the results 
of tests of this code against the presented exact solution. 



2. Equations of relativistic hydrodynamics 

Equations expressing the conservation of the energy and momentum in 
the relativistic hydrodynamics can be written as 

9/.^^" = 0> (1) 

where 



p,i/ 



Tf"" = {p + p)u''v!' +pri' 

is the energy-momentum tensor of the perfect fluid. In this formula p de- 
notes the energy density, p is the pressure, u^ are the components of the 
four-velocity of the fluid, and r/'^^ = diag(— 1, +1, +1, +1) is the metric ten- 
sor of the Minkowski space-timqj. 

For barotropic equations of state of the form p = p{p) the above equa- 
tions constitute a complete set of relevant equations of hydrodynamics. 
However, if the equation of state depends explicitly on the rest-mass density 
n (which is the case for the perfect gas equation of state) we also have to 
take into account the continuity equation, i.e., 

d^inu^ = 0. (2) 

Thus, the set of the equations of hydrodynamics is different for those two 
cases, and this difference becomes especially important if we allow for dis- 
continuous solutions. In the following we will restrict our attention to the 
ultra-relativistic equation of state of the form p = kp. 

Since we aim at solving the general Riemann problem, it is convenient 
to express equations ([T]) in the form 

dtU + diF' = 0, (3) 

where 

V=(^{p + p)W^-p,{p + p)W\\ip + p)W\'^,ip + p)W\^Y , (4) 
and 

F* = (^{p + p)W\\{p + p)W\'v^ + 6'^p, 

{p + p)W^v'v^ + 6'^p, {p + p)W^v'v^ + 6'^p^ ^ . 

Here W denotes the Lorentz factor W = u^, and u* are the components of 
the three- velocity v* = u^/W (note that W = 1/^1 — ViV^). 



^ We will work in Cartesian coordinates in this paper. We will also assume a convention 
in which Greek indices refer to the space-time dimensions /i = 0, 1, 2, 3 while Latin 
ones are reserved for spatial dimensions only i — 1,2, 3. 



dt{{p + p)W^-p) 


+ 


dAi.P + P)W^v^) 


= 0, 


dt{{p + p)W^v^) 


+ 


d^,{{p + p)w^{v-f+p) 


= 0, 


dt{{p + p)w\y) 


+ 


d^{{p + p)W^v-vy) 


= 0, 


dt{{p + p)W\') 


+ 


3.((/3 + p)W^VV) 


= 0. 



3. Solutions of the Riemann problem 

We will now search for the solutions of the Riemann problem for equa- 
tions ([1]) . Without loss of generality, we can assume that the initial discon- 
tinuity is perpendicular to the x axis, and it is located at x = 0. Thus, due 
to the translational symmetry of the initial data in the directions y and z, 
equations ([T|) can be reduced to 



(5) 



The structure of solutions of the Riemann problem for the set of equa- 
tions ([5]) is similar to that obtained in the Newtonian case. In general, the 
initial discontinuity can decay into three kinds of elementary waves (not all 
of them have to be present in the solution) propagating along the x direction 
and separated by some constant states. One of those waves, the so-called 
rarefaction wave, is a smooth self-similar solution of ([5]). The other are two 
discontinuities: a shock wave and a contact discontinuity. The distinction 
between them is based on the behavior of the pressure p and velocity w^', 
which can be discontinuous at a shock wave and have to be equal on both 
sides of a contact discontinuity. In fact, the only quantities that can ex- 
hibit a jump at the contact discontinuity in the case of the ultra-relativistic 
equation of state are v^ and v^ . 

Denoting the Riemann states, that is the initial data for x < and 
X > by L and R respectively, we can symbolically write the solution 
of the Riemann problem as LR — )• L>V^L*Ci?*W^i?. Here the subscript 
arrows refer to the direction from which the particles of the fluid enter the 
wave W, which, in turn, can be either a rarefaction wave 7^ or a shock wave 
S (four different wave patterns corresponding to yV^(^) = 5^(^), 7^^(^) 
are possible). The symbol C denotes a possible contact discontinuity and 
L*, i?* are some intermediate constant states. 

The exact form of the solution can be found by considering the depen- 
dence of the energy density p on the velocity v^ behind the waves W^ and 
W-).. (The exact form of this function depends on the states L and R in 
front of the waves.) Since at the contact discontinuity only v^ and v^ can 
be discontinuous, we must have PL,{vf^^) = PR*{vr,)- It turns out that the 
solution of this equation allows us to establish the solution in the interme- 
diate states L^ and R^, and the precise form of W^(^). In the forthcoming 
sections we will derive the required relations p{v^) for both rarefaction and 
shock waves. 



3.1. Rarefaction wave 

Let us consider a rarefaction wave, that is a smooth solution of (0) 
depending on x and t through ^ = x/t. Under such an assumption the 
equations ([5]) can be reduced to 






ei((p + p)W^V) = fA{p + p)W^v-v-'^ 



(6) 



Non-trivial solutions to the above homogeneous set of ordinary equations 
exist only if the Wronskian of the whole system vanishes, i.e., when ^ are the 
eigenvalues of the Jacobian 9F^/5U. These eigenvalues can be computed 
yielding 

z;^(l -k)± y/kJ{l - ViV^) (1 - v^v^k - {v^Y{l - k)) 

Co = v"", C± = ^- , 

1 — ViV^k 

(7) 

where the eigenvalue ^o is twofold degenerate |i)J. In the special case with 

yy = v^ = 0, the eigenvalues C± can be expressed as 

v^ ± y/k 



1 lb \fkv^ 

where Vk can clearly be identified with the local speed of sound (in the lin- 
earized picture the eigenvalues of 9F^/9U can be interpreted as the speeds 
of propagation of acoustic signals) . 

It can be shown next that for ^ ^ v^ one have 

p'^Wyy = const, p'^Wv'' = const, 

where k = k/{l + k). An assumption that ^ = v^ leads us to the contact dis- 
continuity, which will be discussed later. Introducing the tangential velocity 
v^ = \/{v^)''^ + {v^)'^, we can show that f* = aW~^p~'^, where a denotes a 
constant. 

At this stage equations ^ can be integrated. For a = (no tangential 
velocities) one obtains 

1 -\- V^\ 2 J^ 

= Cip^, 



while for non-zero tangential velocities we get 



1 






^ Vk + ^1 + {I - k)a'^P^' 

where Ci and C2 are integration constants, that can be computed knowing 
the state in front of the wave (for more details concerning this solution see 

3.2. Shock wave 

The solution in the form of a shock wave can be obtained from the 
appropriate Rankine-Hugoniot conditions. For the set of equations ([1]) they 
can be expressed as 

where n'^ is the unit vector normal to the surface of discontinuity. The 
symbol [[•]] is used to denote the difference between the limits of a given 
function at both sides of the discontinuity. As we are primarily interested 
in establishing the state behind the shock wave basing on the known state 
in front of it, we assume a notation in which the value of a given quantity in 
front of the shock wave is denoted with a dash, while an unaltered symbol 
is reserved for the value behind the discontinuity. In this notation a jump 
of a given quantity / is denoted as [[/]] = f — f (a similar convention was 
adopted in jTTj). 

For the discontinuity surface being a plane perpendicular to the x axis, 
the components of the normal vector n^ can be written as n^ = Ws{Vs, 1, 0, 0), 
where Ws = l/\/l — V^. Here Vs has the interpretation of the shock wave 
velocity. 

The Rankine-Hugoniot conditions can be now written as 

[[pW^-Kp]]Vs = [[pW^v^]], 

[[pWV]]y, = [[pWHvn^ + Kp]] , 

[[pw\y]] Vs = [[pw^v^'vy]] , ^^ 

[[pW'^v']] Vs = [[pW^v^v']] . 

For v^ = v^ = the only physical solution of the above algebraic set of 
equations is given by 



pl + G + V(l + 0)2-1 



where 6 = W'^W^{v'' 
be expressed as 



r,x\2 



/(2k(1 — k)). The shock wave velocity Vg can 



V, 



pW^v"" I pW^ - tip 



For non-zero tangential velocities v* the speed of propagation of the 
shock wave can be obtained as the root of the following cubic equation 



:i - ^"K 



(1 - v^Vs){l - v^Vs) - \{v^ - VsW 

-t\2 



Vs 



{vy{l-v^Vs){l-Vf) 



0. 



The value of Vs can always be computed with the help of Cardano's formu- 
lae, but a Newton-Raphson scheme may be more efficient in numerical ap- 
plications. The above equation was derived from ([8]) under the assumption 
that Vs ^ v^. The case with Vs = v^ , corresponds to a contact discontinuity. 
Having obtained a solution for Vg we can express the value of the energy 
density behind the shock wave as 

^ {v--Vs){l-v-Vs){l-v-Vs) ■ ^ ^ 

The solution is completed by the expression for the tangential velocity 

{Vs - v'^fp^W^v'f 



,,t^2 



P^W\Vs 



,x\2 



3.3. Contact discontinuity 

Apart form the shock waves described above, equations ([8]) admit non- 
trivial solutions for which Vs = v^, that is, the discontinuity is moving with 
the flow of the fluid. In this case we have v^ = v^{= Vg) and p = p. There 
is, however, no constraint on the tangential velocity, and it can exhibit an 
arbitrary jump. Such a discontinuity is called a contact one. 

W should also note that in this respect the hydrodynamics with the 
ultra-relativistic equation of state differs qualitatively from that of the per- 
fect gas equation of state. In the latter case a contact discontinuity can also 
be present due to the jump in the rest mass density and the specific internal 
energy. 

3.4. Solutions of the Riemann problem 

Whether we are dealing with the rarefaction or a shock wave can be 
distinguished basing on the ratio of the pressure p in front of the wave to 
the pressure p behind the wave. For p > pwe have a rarefaction wave, while 



the condition p < p is characteristic for a shock wave [12j. It follows from 
the results presented in the preceding sections that this distinction can be 
translated into a condition concerning the velocity v^. 

Let p = 5_>(^)(f^) be the energy density behind the shock wave under- 
stood as a function of the post-shock velocity v^ and given by equation ([9]) 
(here again the arrows refer to the direction from which the fluid enters the 
wave) . Let us also introduce a similar function for the rarefaction wave and 
denote it by p = 7^_^(^)(f^). The expressions for the energy density behind 
an arbitrary wave W_5.(^) can be now written as 



for a right moving wave, and 
p = W^iv^ 






for a left moving one. The symbol v^ used here denotes the normal velocity 
in front of a wave. In order to show the influence of the tangential velocities 
on the solution, we depicted several of such functions on Fig. [H plotting the 
graphs for different values of u* for both right and left-moving waves. 

The strategy of finding the solution of the Riemann problem can be now 
described as follows. We consider a function p = W<_(f^) for the state L 
and p = W->(w^) for the state R, where the states L and R represent the 
data in front of the waves. The intersection of the graphs of those functions 
gives the values of v^ and p^, common for both intermediate states L* 
and i?*, and also identifies the character of both waves (the so-called wave 
pattern) . To complete the solution one only needs to establish the velocities 
with which the fronts of the waves, and the tail of the rarefaction wave (if 
present) propagate. The shock wave moves with the velocity Vg, which can 
be easily computed once the value of v^ has been obtained. The velocity 
of the head of the rarefaction wave is given by (,± (plus for TZ^ , minus for 
7?.^) calculated for the suitable Riemann state. The location of the tail of 
the rarefaction wave can be computed from the condition that the energy 
density in the wave should reach the value of p* . 

A sample solution obtained in this way is shown on Figs. [2] and [3j As 
the initial states for this example we chose pL = 1, f£ = 1/2, v\^ = 1/3 
for the left state and pR = 20, i^^ = 1/2, u^ = 1/2 for the right one. The 
equation of state was assumed to be p = p/3. 

4. Description of the numerical code 

We will now describe a numerical code that has been implemented in 
order solve equations ([3]) in (3 -|- 1) dimensions. 
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Fig. 1. The dependence of the energy density p on the velocity v^ behind the wave 
for the uhra-relativistic equation of state with k = 1/3. Different curves refer to 
values of the tangential velocity w* in front of the wave equal to 0, 0.5, 0.8, and 
0.865. The velocity v^ in front of the wave is equal 0.5, and the density p was set to 
1. Increasing curves correspond to the right moving waves, while decreasing ones 
to the left moving waves. 



The construction of the code is similar to the one described in [13] . The 
vector of conserved quantities U given by Q is discretized on a Cartesian 
grid of cells and its time derivative is computed according to the following 
method of lines 



dt 



d'^i,j,k i_ fj^x _ F^ ^ f-py _ T^y 

~ Ax V *+V2,i,fc ' i-l/2,j,kj ^y [^ i,j+l/2,k ^ i,j~l/2,k 






Here indices i, j, k number the cells of the grid in the directions x, y and 
z, while F|j_^y2,j,fc' ^l,j±i/2,k ^^^ '^i,j,k±i/2 ^^^ numerical fluxes at cells' 
interfaces. With this discretization the values of U,- 



'ij^k are evolved with the 
standard Runge-Kutta method of the second order. The time step for this 
evolution is limited by the Courant criterion, in which it is convenient to 
assume the speed of light c = 1 as the upper bound for the velocity of any 
signals. 

The numerical fluxes are computed following a method developed by 
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Fig. 2. The solution of the Riemann problem for t — 1. The left initial state is 
given by p^ = 1, w£ = 1/2, v^ = 1/3 and the right state by pB. — 20, vf^ — 1/2, 
v'h - 1/2. 



Harten, Lax and van Leer [14J (HLL method, hereafter) and proposed for 
the relativistic hydrodynamics by Schneider et al. [15]. A flux between two 
cells characterized by the states Ul and \Jr is given by 



CmaxF(UL) + CmmF(Uij) - CmaxCmin (U R - V l) 



+ Cn 



.L,R 



where Cmax = max{0, Aq'", A^' } and Cmin = — min{0,Ao' , A^' }. Here 
Aq' and A_|_' denote the eigenvalues of the Jacobians dY^ /d\J as given by 
the formula ([7]) and computed for the states Ul and U/j respectively. 

The so-called minmod reconstruction method of Van Leer [16] was used 
to improve the spatial accuracy of the scheme. In this approach the states 
^L,R used to compute numerical fluxes are obtained as follows. We first 
define the "minmod" function 





r a 


if \a\ < \h\, 


ah > 0, 


minmod(a, h) = < 


^ 


if \a\ > \h\, 


ab > 0, 







if a6 < 
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Fig. 3. Time snapshot of the sohition of the Riemann problem with the same initial 
data as on Fig. [21 The solid line corresponds to the velocity w*, while the dotted 
one depicts w^. 



U,+i-U, Ui-U,_i 



and the slope limiters 

S,; = minmod 

Then the expressions for left and right states used to compute the numerical 
flux at the interface i + 1/2 are given by 

U; 



2 — Uj + Sj ( Xj4_i/2 



Uj+l/2 = Uj+i + Sj+i [xi^i/2 - Xi+i 

The above formulae actually refer to the reconstruction procedure in x di- 
rection. Formulae for other directions are analogous. 

The recovery of primitive variables, i.e., obtaining the values of p, p, v^, 
v'^ and v^ from U is straightforward, and, unlike in simulations with the 
perfect gas equation of state, it does not involve a numerical solution to 
a nonlinear algebraic equation. In the latter case this is usually obtained 
using an iterative Newton-Raphson scheme. This fact results in a robust 
performance of the code for the ultra-relativistic equation of state, as com- 
pared to the standard case. We omit the exact formulae for the recovery of 
the primitive quantities, as they follow directly from ([H). 



12 



22 

20 

18 

16 

14 

12 

10 

8 

6 

4 



exact solution 
numerical solution 




-0.5 





X 



0.5 



Fig. 4. Sample numerical solution for the energy density (points) plotted over the 
exact solution from Fig. [5] (solid line). The numerical solution was obtained with 
the spatial resolution of 800 zones per unit length. 



Boundary conditions are implemented by adding to the numerical do- 
main 2 "ghost" zones in each direction, following the standard procedure 
used in numerical hydrodynamics (see e.g. |I7|). The method of "ghost" 
zones is also used to split the grid into parts in the parallel Message Passing 
Interface (MPI) implementation of the code. In this case the values of the 
"ghost" zone variables are obtained from adjacent parts of the grid. 



5. Sample code tests 

We used the exact solution presented on Figs. [2] and [3] to test the code. 
To this end we performed a couple of 2-dimensional runs, so that the tangen- 
tial velocity was also evolved, reaching the time t = 1. We used equidistant 
grids of 100, 200, 400, 800 and 1600 zones per unit length. In each case 
the Courant factor was set to 0.1. Results of a run with the resolution of 
800 zones per unit length are shown on Figs. U] and [5j The data depicted 
on these graphs come from one y = const slice of the grid. We see that 
the resolution of discontinuities is worse for the tangential velocities, most 
probably due to non-negligible numerical viscosity of the scheme. 
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Fig. 5. Sample numerical solutions for the velocities u^ and u* (points) plotted 
over the exact curves from Fig. [31 The numerical solutions were obtained with the 
spatial resolution of 800 zones per unit length. 



Zones / unit length 



P 



V 



yV 



100 
200 
400 
800 
1600 



3.1 
1.7 
9.2 

4.7 
2.5 



"10^ 
10-1 
10-2 
10^2 

10-2 



1.6 
9.1 
6.2 

2.8 
1.7 



10-3 
10-3 
10-3 

10-3 



1.8 
1.0 
6.6 

4.1 
2.5 



10-2 

10-3 
10-3 
10-3 



Table 1. Estimates of absolute errors of sample numerical runs for the Riemann 
problem illustrated on Figs. [H and |31 The errors were computed in the Li norm on 
the interval x G [—1,1] with respect to the exact solution. 



Table [T] gives the estimates of the absolute numerical errors computed 
with respect to the Li norm on the interval x G [—1,1]. The exact solution 
was taken as the reference for the calculation of these errors. The linear 
convergence in ah variables is characteristic for test problem involving strong 
discontinuities. 

We should also point out that the HLL Riemann solver implemented 
in the scheme is an approximate one and it does not exploit the analytical 
solution presented in this chapter. Thus the presented test and the imple- 
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mentation of the code are independent — we are not testing a given solution 
against itself. 

6. Summary 

We have presented the exact solution of the Riemann problem in the rela- 
tivistic hydrodynamics with the ultra-relativistic equation of state, in which 
the fluid is allowed to move with arbitrary velocities, also tangent to the sur- 
face of the initial discontinuity. We have also described a 3-dimensional nu- 
merical scheme for solving relativistic equations of hydrodynamics with the 
ultra-relativistic equation of state. Tests of the implemented code against 
the obtained analytic solution show satisfactory convergence and accuracy 
allowing to hope for future applications of the code. 

We are also convinced that the analytic solution discussed in this article 
can also be used to construct an exact Riemann solver for the 3-dimensional 
hydrodynamic codes with the ultra-relativistic equation of state. 

We also believe that this solution can find its applications outside strictly 
numerical hydrodynamics. It was already mentioned in the beginning that a 
similar solution for the perfect gas equation of state can be used to explain a 
boosting mechanism present in relativistic jets [TO]. The required "boosting" 
property is also exhibited by solutions of this paper. This fact can be 
noticed, for instance, by looking at the solution from Fig. El where the 
tangential velocity in the rarefaction wave reaches values much greater than 
any of velocities present in the initial states. 
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